Ancient DNA and osteological analyses of a unique paleo-archive reveal Early Holocene faunal expansion into the Scandinavian Arctic

Paleo-archives are essential for our understanding of species responses to climate warming, yet such archives are extremely rare in the Arctic. Here, we combine morphological analyses and bulk-bone metabarcoding to investigate a unique chronology of bone deposits sealed in the high-latitude Storsteinhola cave system (68°50′ N 16°22′ E) in Norway. This deposit dates to a period of climate warming from the end of the Late Glacial [~13 thousand calibrated years before the present (ka cal B.P.)] to the Holocene thermal maximum (~5.6 ka cal B.P.). Paleogenetic analyses allow us to exploit the 1000s of morphologically unidentifiable bone fragments resulting in a high-resolution sequence with 40 different taxa, including species not previously found here. Our record reveals borealization in both the marine and terrestrial environments above the Arctic Circle as a naturally recurring phenomenon in past periods of warming, providing fundamental insights into the ecosystem-wide responses that are ongoing today.


DNA extraction and amplification
Pre-digestion and DNA extraction protocols followed Lord et al. (100).Samples were pre-digested in 700 ųl extraction buffer (Urea 1M, EDTA 0.5M at pH8) with 15 ųl proteinase K (18 mg/ml) at 55 o C for 30 minutes.Samples were subsequently digested overnight on a nutator at 55 o C with fresh extraction buffer and proteinase K. Following digestion, 500 ųl of supernatant was concentrated to ~100 ųl using an Amicon Ultra-0.5 Centrifugal Filter Unit (30kDa MWCO), purified with a MinElute PCR Purification Kit (QIAGEN) and eluted in 100 ųl elution buffer (EB buffer).Primers targeted a region of the 16S rRNA gene of the mitochondria for amplification in Mammals and Fish, and a 12S region for birds.Mammalian DNA was amplified using the Mamp007 primers [5'-CGAGAAGACCCTATGGAGCT-3', 5'-CCGAGGTCRCCCCAACC-3'; (102)], fish DNA was amplified using the Fish16S primers (5'-TACCAAAAACATCGCCTCYTG-3', 5'-CATTTAAAAGACAAGTGATTRCG-3'; this study, designed by Laura Epp) and bird DNA with the Aves12S primers [5'-GATTAGATACCCCACTATGC-3', 5'-GTTTTAAGCGTTTGTGCTCG-3'( 101)].Samples with unknown bone fragments were amplified with all three primer pairs.PCR reactions (25 ųl) contained 2 units of AmpliTaq Gold DNA Polymerase (Applied Biosystems), 15 mM (Tris-HCL), 150 mM (KCl), 10X amplification buffer, 0.2 mM of each dNTP, 2.5 mM MgCl2, 0.8 ug BSA, 0.2 ųM of each primer and 5 ul of template DNA.For amplifications with Mamp007, two units of blocking primer MamP007_B_Hum1, 5′-GGAGCTTTAATTTATTAATGCAAACAGTACCC-3′ were added to reduce the amplification of human DNA (102,105).Forward and reverse primers were tagged with a unique 8 or 9 bp barcode at the 5' end to allow for multiplexing, with each primer pair having the same tag.The list of tags used can be found in Taberlet et al. (106).PCR cycling conditions were as follows: 10 minutes at 95 o C, 40 cycles at 95 o C for 30s, annealing temperature for 30s, 72 o C for 1 minute, and a final extension of 10 minutes at 72 o C. Annealing temperature was 55 o C for the mammal and bird primers, and 50 o C for the fish primers.Negative controls were included for all DNA extractions and PCRs.All laboratory protocols up to PCRs were carried out in the dedicated ancient DNA laboratory at the University of Oslo following standard protocols to minimize contamination.Amplified products were visualized on agarose gels (1.5%) and pooled according to the strength of the bands (three categories were recognized: strong, medium and weak).Pools were purified using the MinElute PCR Purification Kit (QIAGEN) and concentrations of pools were measured using the Qubit 2.0 dsDNA BR Assay Kit (Thermo Fisher).Pools were subsequently combined equimolarly, and libraries were built using the TruSeq DNA Nano library preparation kit and sequenced on the Illumina NovaSeq 6000 S4 (150 bp PE) at the Norwegian Sequencing Center.

Data
was processed using the OBITools package v.1.2.12 [https://pythonhosted.org/OBITools/index.html; (99)] following https://pythonhosted.org/OBITools/wolves.html.Forward and reverse reads were assembled using illuminapairedend and samples were assigned with ngsfilter.Subsequently, reads were discarded if their quality score was < 40, tags had < 100% match, primers had >3 mismatches, or read length was < 20bp.Unique reads were merged and obiclean was used to remove singletons and identify likely PCR and/or sequencing artefacts by applying a 5% threshold ratio to reclassify "internal" sequences to their corresponding "head".For taxonomic assignment sequences were compared to reference libraries using ecotag.Reference libraries were built for each of the different primer pairs (Mamp007, Fish16S and Aves12S) by performing an in-silico PCR with ecoPCR (104)  Following analyses with ecotag sequences were filtered in R v.4.3.0 (https://www.r-project.org/) to minimize misidentification by removing sequences with < 95% identity for fish and birds and < 98% for mammals, taxa with a total read count under 200 and PCR replicates with less than 100 reads total.No reads were assigned to the common contaminants Bos sp. or Gallus sp., and following filtering no reads assigned to the genera Homo, Ovis or Sus remained in the dataset.Postfiltering, one extraction negative control had reads identified to the Leuciscinae fish family.One sample contained reads identified to Leuciscinae (261 reads), which were conservatively removed.Post-filtering no sequences remained in any of the other negative controls.
For comparison of replicate PCRs only samples where both PCR replicates contained reads after filtering and analyses were considered.This led to a total dataset of 22 samples for the Mamp007 primer, 4 samples for the Aves12S primer, and 35 samples for Fish16S.A presence/absence matrix of each PCR repeat for a given sample was generated.The number of extra taxa was calculated, and a mean and standard deviation generated per taxonomic group (mammals, birds, fish).Analyzing two PCR replicates allowed for the identification of an additional 0.49 ± 1.07 taxa in Aves12S, 1.68 ± 2.61 in Mamp007, and 1.75 ± 1.71 taxa in Fish16S.Considering the similarity between replicates we merged taxonomic identifications from both PCR replicates for further downstream analyses (Table S9).
We also analyzed the effect of extracting DNA from multiple subsamples from a single bulk-bone sample, specifically assessing the benefit of two or three subsamples instead of just one.A total dataset of 13 samples were used for analyzing the effect two subsamples versus one subsample, and 10 samples were used for analyzing the effect of analyzing three subsamples versus two subsamples as well three subsamples versus one subsample.To avoid any impact of PCR replicates into this test, only PCR1 was used in this analysis.We generated a presence/absence matrix as described above for the analyses of PCR replicates and calculated the average number of additional taxa that was detected when adding subsamples (Table S10).Analyzing two subsamples versus one identified on average 2.92 ± 3.59 more taxa.Analyzing an additional sample (i.e., three subsamples versus two subsamples) added on average 2.70 ± 3.16 more taxa (Table S10).Finally, analyzing three subsamples versus just one lead to the identification of 5.00 ± 4.24 more taxa.We conclude that adding subsamples leads to the identification of more taxa and may be beneficial when sufficient material is present.For the current study we merged identifications from different subsamples per layer.

Adjustments of taxonomic identification
The complete taxa list was reviewed by several taxonomic experts with a good understanding of Northern European fauna.Anne Karin Hufthammer, professor at the University Museum of Bergen, Liselotte M. Takken-Beijersbergen, senior engineer at the University Museum of Bergen, Thijs van Kolfschoten, professor emeritus in archeology at the University of Leiden, reviewed the mammals.Ingvar Byrkjedal, professor emeritus and Nicolas Straube, associate professor and curator of ichthyology from the Department of Natural History at the University Museum of Bergen reviewed the fish.Samuel J. Walker at the University of Oslo reviewed the birds.Sequences were thereafter assessed with NCBI Nucleotide BLAST tool as a complementary analysis for taxa that were considered unlikely or unexpected to be present.When multiple sequences were identified to the same taxon, the sequence with the highest identity and highest abundance was selected for analyses with BLAST.A subset of taxonomic identifications was adjusted based on the assessment of the taxonomists and the BLAST results (Table S6).After any manual taxonomic adjustments, the taxa lists from aDNA BMM and osteological analyses were combined.S1.Modelled using OxCal v.4.4.4 (93).DRY6, in blue, gave a terrestrial reading and was therefore modelled using the IntCal20 calibration curve for the Northern Hemisphere (96).All other dates were modelled on the Marine20 calibration curve (94) with a ΔR offset of -100 ± 37 (95).
Table S1.Radiocarbon dates for Nygrotta.Uncalibrated radiocarbon ages and calibrated ages at 95.4% probability.δ 13 C values are given to justify marine or terrestrial models used to calibrate.BP = before present.Modelled using OxCal [Fig.S1; (93)], DRY6 gave a terrestrial reading and was therefore modelled using the IntCal20 calibration curve for the Northern Hemisphere (96).All other dates were modelled on the Marine20 calibration curve (94) with a ΔR offset of -100 ± 37 (95).

Table S5. Number of reads per taxa identified by bulk-bone metabarcoding presented by stratigraphic layer (A, B, C).
Original taxonomic identifications are presented, and identifications that were later manually adjusted or removed are marked with an asterisk (see Table S6).Table S6.Manually adjusted or removed taxonomic identifications.Justifications for manual adjustment and removal of taxa identified from the bulk-bone metabarcoding data using the OBITools package v. 1.2.12 (103).

EcoTag ID Adjusted ID Justification
Centrocercus sp.Phasianidae Distributed in North America (108).We adjusted this identification to the family level of Phasianidae.

Centrocercus urophasianus
Phasianidae Distributed in North America (108).We adjusted this identification to the family level of Phasianidae.
Dendragapus sp.Phasianidae Distributed in North America (108).We adjusted this identification to the family level of Phasianidae.
Tympanuchus sp.Phasianidae Distributed in North America (108).We adjusted this identification to the family level of Phasianidae.

Lepus sp.
Lepus timidus Results of the NCBI-Blast gave a 100% match with five species of Lepus.As Lepus timidus (mountain hare) is the only species present in Fennoscandia, we shift this identification to species level.
Results of the NCBI-Blast gave a 100% match on both query cover and sequence identity with wildcat Felis silvestris (European wildcat), African wildcat (Felis lybica) and Felis catus (domestic cat), and no match with the native Lynx lynx (Eurasian Lynx) despite it being in the NCBI database.We adjust the identification to Felis sp.
Ipnops sp.Removed Ipnops sp. are deep sea fish and do not occur in the Norwegian Sea.The family Ipnopidae (deep-sea tripod fishes) are also not present in the Norwegian sea or surrounding waters.NCBI-Blast showed no 100% query cover matches.The sequence was amplified by primer pair Fish16S, does not match any common contaminants and its identification may instead reflect incomplete coverage of the reference database.Adjusting the taxonomic identification would lead to an identification beyond the family level.We removed this identification.

Esox masquinongy Esox lucius
Native to North America (109).Results of NCBI-Blast were inconclusive, with no 100% query cover matches.However, as Esox lucius (northern pike) is the only Esox species in Fennoscandia we adjust the identification to this species.111))] we adjusted this identification to the species L. limanda.

Sebastiscus marmoratus
Sebastinae Tropical species distributed in the Western Pacific.Interestingly one has been found in Norwegian waters (112) and this seems to be the only recording of the species in the North Atlantic.How it reached the North Atlantic is unknown but likely as larva or fry via a ballast water.Nevertheless, we adjusted the identification to the family of Sebastinae which has many taxa in the Atlantic.

Artedius lateralis
Cottidae Native to the North Pacific (113).NCBI-Blast resulted in 13 matches at 100% query cover and 98.75% sequence identity, 10 of which are distributed in the North and East Pacific.Three are present in the Norwegian Sea: Myoxocephalus Scorpius (shorthorn sculpin), Myoxocephalus quadricornis (fourhorn sculpin) and Gymnocanthus tricuspis (Arctic staghorn sculpin).We adjusted this identification to the family level Cottidae.

Enophrys diceraus
Cottidae Native to the North Pacific (113).NCBI-Blast results gave a 100% query cover and 98.75% sequence identity match with 16 species, the majority of which are distributed in the North and Northwest Pacific.Three species are currently found in the Norwegian Sea: Gymnocanthus tricuspis (Arctic staghorn sculpin), Myoxocephalus Scorpius and Myoxocephalus quadricornis.We adjusted this identification to the family level Cottidae.
Native to the North Pacific (113).NCBI-Blast results show a 100% query cover and 98.67% sequence identity match with two species native to Norway: Anarhichas lupus (Atlantic wolffish) and Anarhichas minor (spotted wolffish).Anarrhichthys ocellatus (wolf-eel) was also a 98.67% match, however it only occurs in the North Pacific (113).We adjusted this identification to the genus Anarhichas.

Microhylinae Removed
Microhylinae are a subfamily of Mycrohylidae (narrow-mouthed frogs), few of which occur in the Northern Hemisphere and none of which occur in Europe.The sequence was amplified by primer pair Fish16S, it does not match any common contaminants and its identification may instead reflect incomplete coverage of the reference database.Adjusting the taxonomic identification would lead to an identification beyond the family level.This identification was removed.
Nasikabatrachus sp.Removed Endemic to India and the only genus in its family.NCBI Blast gave a match with 100% query cover and 97.26% sequence identity with two species endemic to the African continent (Amnirana galamensis and Amnirana albolabris).The sequence was amplified by primer pair Fish16S , it does not match any common contaminants and its identification may instead reflect incomplete coverage of the reference database.Adjusting the taxonomic identification would lead to an identification beyond the family level.This identification was removed.
Table S7.Taxa identified by BBM and osteology by stratigraphic and mechanical layer for layer B (B1, B2, B3).Layer A is dated to 5948-5611 cal BP, layers B1, B2 and B3 are dated to 9534-9006 cal BP and layer C is dated to 13,069 -12,744 cal BP.
Table S8.Sample information of the 20 BBM samples analyzed.Layer presents the stratigraphic (letter) and mechanical (number) layer that the sample comes from.The number of subsamples from which DNA was extracted, the total weight of the available material and the total number of bone fragments in the total sample are given.DNA was extracted from up to three subsamples of about 110 mg where that amount was available from a total of 47 subsamples.

Figure S1 .
Figure S1.Modelled plot of radiocarbon dates for Nygrotta.Information on each sample, including the uncalibrated and calibrated dates are given in TableS1.Modelled using OxCal v.4.4.4(93).DRY6, in blue, gave a terrestrial reading and was therefore modelled using the IntCal20 calibration curve for the Northern Hemisphere(96).All other dates were modelled on the Marine20 calibration curve (94) with a ΔR offset of -100 ± 37(95).

Table S2 . Bone element distribution by taxonomic group as identified by osteology.
All fragments identified to bone element are included from all layers combined.

Table S3 . Number of Identified Specimens (NISP) from Nygrotta per layer.
16All taxa identified by osteology and their corresponding NISP figure, which was calculated taking into account all bone elements (including vertebrae and phalanges).Unidentified fragments were identified to class (Aves, Mammalia and Pisces) where possible and into unidentified if class could not be determined.Ages are given in calibrated years before present (cal BP).

Table S5 continued. Number of reads per taxa identified by bulk-bone metabarcoding presented by stratigraphic layer
Original taxonomic identifications are presented, and identifications that were later manually adjusted or removed are marked with an asterisk (see TableS6).

Table S6 continued. Manually adjusted or removed taxonomic identifications.
(103)fications for manual adjustment and removal of taxa identified from the bulk-bone metabarcoding data using the OBITools package v. 1.2.12Boyer et al.Boyer, Mercier, Bonin, Le Bras, Taberlet and Coissac(103).

Table S9 . Average number of extra taxa detected by doing a second PCR replicate, per PCR primer.
Only replicate pairs with reads remaining after filtering were considered.22 samples were used for this analysis for the Mamp007 primer, 4 for Aves12S, and 35 samples for Fish16S.

Table S10 . The average number of extra taxa detected when extracting DNA from two or three subsamples
. Only PCR1 was used for this test and only subsamples with reads postfiltering were included.